The assignment is to use numerical methods to generate pressure and friction force predictions for the laminar flow around three geometries:
Specifically, for the circular cylinder you will:
For the jukowski foil you will:
The mystery geometry will require a similar computation. Correct results for this geometry will demonstrate independent excellence and will be marked accordingly.
For the circular cylinder you will:
Firstly, I use the numpy to caculate the velocity u and v: $$u(x,y) = \frac{\gamma}{2\pi}\left[\tan^{-1}\left(\frac{x-S}y\right)-\tan^{-1}\left(\frac{x+S}y\right)\right].$$
$$v(x,y) =\frac{\gamma}{4\pi} \log\left(\frac{(x+S)^2+y^2}{(x-S)^2+y^2}\right).$$It has the code:
In [1]:
import numpy
from matplotlib import pyplot
%matplotlib inline
#Import the required functions from VortexPanel.py and BoundaryLayer.py
from VortexPanel import Panel, solve_gamma, plot_flow, make_circle
In [2]:
pyplot.figure(figsize=(10,6))
def c_p(gamma): return 1-gamma**2
def C_P(theta): return 1-4*(numpy.sin(theta))**2
N_resolution=[32,64,128]
for m in range(3):
circle = make_circle(N_resolution[m])
solve_gamma(circle)
coeff = numpy.zeros(N_resolution[m])
for i,p_i in enumerate(circle): coeff[i]=c_p(p_i.gamma)
pyplot.plot(numpy.linspace(1,N_resolution[m],N_resolution[m])/N_resolution[m]*numpy.pi*2,coeff)
pyplot.legend(["resolutions(32)","resolutions(64)","resolutions(128)"])
pyplot.title('$C_p$ of the circular cylinder under different resolution')
pyplot.xlabel('phase angle /rad')
pyplot.ylabel('$C_p$')
pyplot.plot(numpy.linspace(0,numpy.pi*2,180),1-4*numpy.sin(numpy.linspace(0,numpy.pi*2,180))**2,lw=2, c='k')
Out[2]:
The black line is the exact potential flow solution, the figure shows when we increase the resolutions which means use more numbers of pannels to model can make the result more accurate.
(2)Compute the friction drag coefficient $C_F$, and an estimate for the pressure drag coefficient $C_D$ for the same three resolutions at $Re_D=10^5$.
In [3]:
def pohlF(eta):
return 2*eta-2*eta**3+eta**4
def pohlG(eta):
return eta/6*(1-eta)**3
def pohl(eta,lam):
return pohlF(eta)+lam*pohlG(eta)
In [4]:
def disp_ratio(lam): return 3./10.-lam/120.
def mom_ratio(lam): return 37./315.-lam/945.-lam**2/9072.
def df_0(lam): return 2+lam/6.
In [5]:
def g_1(lam): return df_0(lam)-lam*(disp_ratio(lam)+2*mom_ratio(lam))
In [6]:
from scipy.optimize import bisect
lam0 = bisect(g_1,-12,12) # use bisect method to find root between -12...12
print 'lambda_0 = ',lam0
In [7]:
def ddx_delta(Re_d,lam):
if Re_d==0: return 0 # Stagnation point condition
return g_1(lam)/mom_ratio(lam)/Re_d # delta'
In [8]:
def heun(g,psi_i,i,dx,*args):
g_i = g(psi_i,i,*args) # integrand at i
tilde_psi = psi_i+g_i*dx # predicted estimate at i+1
g_i_1 = g(tilde_psi,i+1,*args) # integrand at i+1
return psi_i+0.5*(g_i+g_i_1)*dx # corrected estimate
In [9]:
def g_pohl(delta_i,i,u_e,du_e,nu):
Re_d = delta_i*u_e[i]/nu # compute local Reynolds number
lam = delta_i**2*du_e[i]/nu # compute local lambda
return ddx_delta(Re_d,lam) # get derivative
In [10]:
def march(x,u_e,du_e,nu):
delta0 = numpy.sqrt(lam0*nu/du_e[0]) # set delta0
delta = numpy.full_like(x,delta0) # delta array
lam = numpy.full_like(x,lam0) # lambda array
for i in range(len(x)-1): # march!
delta[i+1] = heun(g_pohl,delta[i],i,x[i+1]-x[i], # integrate BL using...
u_e,du_e,nu) # additional arguments
lam[i+1] = delta[i+1]**2*du_e[i+1]/nu # compute lambda
if abs(lam[i+1])>12: break # check stop condition
return delta,lam,i # return with separation index
In [11]:
nu=1e-5 # viscosity
N = 32 # number of steps
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta,lam,iSep = march(s,u_e,du_e,nu) # solve!
In [12]:
def half_c_f(lam, nu, delta, u_e):
Re_d = delta*u_e/nu
return df_0(lam)/Re_d
In [13]:
def taw_w(lam, nu, delta, u_e):
if u_e == 0:
return 0
else:
return half_c_f(lam, nu, delta, u_e)*u_e**2
In [14]:
def C_F_calc(tau, sx, N):
#return numpy.sum(tau[:iSep+1]*sx*numpy.pi/(N-1))
return numpy.trapz(tau[:iSep+1]*sx, dx = numpy.pi/(N-1))
In [15]:
for N in N_resolution:
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta, lam, iSep = march(s,u_e,du_e,nu) # solve!
taw = numpy.full_like(delta, 0)
taw = [taw_w(lam[i], nu, delta[i], u_e[i]) for i in range(N)]
sx = numpy.sin(s[0:iSep+1])
print ('When N = ' + '%i' %N)
C_F_circle = 2*C_F_calc(taw, sx, N)/numpy.pi
print ('Circle frictional coefficient = ' + '%0.2e' %C_F_circle)
C_F_flat = 1.33 * numpy.sqrt(nu/numpy.pi)
print("Flate plate: "+'%0.2e' %C_F_flat)
By using the code numpy.trapz to integrate and get the drag coefficient in three resolutions are:
| Resolution | Circle | Flat plat |
|---|---|---|
| 32 | 4.04e-03 | 2.37e-3 |
| 64 | 4.06e-03 | 2.37e-3 |
| 128 | 4.07e-03 | 2.37e-3 |
In [16]:
s_x = numpy.sin(s)
pyplot.figure(figsize=(8,6))
pyplot.plot(s,taw*s_x)
pyplot.scatter(s[iSep+1], taw[iSep+1]*s_x[iSep+1],s=50,c="r")
pyplot.xlabel('$s$',size=20)
pyplot.ylabel(r'$\tau_w s_x$', size=20)
Out[16]:
Estimate for the pressure drag coefficient CD
In [17]:
def cylinder_separation(s,iSep):
c1=C_P(s[0:iSep+1])
c2=C_P(s[iSep])*numpy.ones(s.size-iSep-1)
return numpy.concatenate((c1,c2),axis=0)
In [18]:
N=64
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta, lam, iSep = march(s,u_e,du_e,nu) # solve!
C_P_s= cylinder_separation(s,iSep)
In [19]:
pyplot.figure(figsize=(10,6))
pyplot.plot(s,cylinder_separation(s,iSep))
pyplot.xlabel('AoA alpha (rad)',size=20)
pyplot.ylabel('$Pressure coefficient C_p$', size=20)
pyplot.scatter(s[iSep+1],C_P_s[iSep+1],s=50,c="r")
pyplot.legend(["$C_p$"])
Out[19]:
In [20]:
def coefficient_Cd(N,sy,C_P_s):
return numpy.trapz(C_P_s*sy, dx=numpy.pi/(N-1))
In [21]:
for N in N_resolution:
s = numpy.linspace(0,numpy.pi,N) # distance goes from 0..pi
s_y=numpy.cos(s)
u_e = 2.*numpy.sin(s) # velocity
du_e = 2.*numpy.cos(s) # gradient
delta, lam, iSep = march(s,u_e,du_e,nu) # solve!
C_P_s=cylinder_separation(s,iSep)
print ("The drag coefficient is:")
print coefficient_Cd(N,s_y,C_P_s)
For the jukowski foil you will:
In [22]:
from VortexPanel import Panel,solve_gamma_kutta,plot_flow,make_jukowski,make_circle
N = 64
foil = make_jukowski(N)
In [23]:
#calculate C_L in a range of angle of attack
def C_L_AA(foil,alpha_max,M):
C_L=numpy.zeros(M)
for i in range(M):
alpha=alpha_max*i*numpy.pi/180./(M-1)
solve_gamma_kutta(foil,alpha) # solve for gamma
C_L[i]=lift(foil) # print the lift
return C_L
In [24]:
def CLA(rr,alpha):
return 2.*numpy.pi*(1+4./numpy.sqrt(3)/3.*rr)*numpy.sin(alpha)
In [25]:
def lift(panels):
c = panels[0].x[0]-panels[len(panels)/2].x[0] # length scale
return -4./c*numpy.sum([p.gamma*p.S for p in panels])
In [26]:
pyplot.figure(figsize=(10,6))
N=[32,64,128]
M=101
alpha_max=10.
for p in range(3):
Na=N[p]
foil1 = make_jukowski(Na,dx=0.15,dy=0,dr=0) # make foil
coff1=C_L_AA(foil1,10,M)
pyplot.plot(numpy.linspace(1,M,M)/M*alpha_max,coff1)
M=101
rr=0.15
coff2=numpy.zeros(M)
for q in range(M):
alpha=alpha_max*q*numpy.pi/180./(M-1)
coff2[q]=CLA(rr,alpha)
pyplot.plot(numpy.linspace(1,M,M)/M*alpha_max,coff2)
pyplot.xlabel('Angle of attack /deg')
pyplot.ylabel('$C_L$')
pyplot.legend(["N=32","N=64","N=128","Analytical"])
Out[26]:
As we ca see in the picture: when the resolution gets higher, the lift coefficient which we calculate is closer to the analytical value.
(2) Compute the separation point:
In [27]:
from BoundaryLayer import Pohlhausen, march
def solve_plot_boundary_layers(panels,alpha=0,nu=1e-5):
# split the panels
top_panels,bottom_panels = split_panels(panels)
# Set up and solve the top boundary layer
top = Pohlhausen(top_panels,nu)
top.march()
# Set up and solve the bottom boundary layer
bottom = Pohlhausen(bottom_panels,nu)
bottom.march()
return top,bottom
In [28]:
def predict_jukowski_separation(t_c,alpha=0,N=128):
#Function from separation prediction and altered to return a value.
# set dx to gets the correct t/c
foil = make_jukowski(N,dx=t_c-0.019) #t_c-0.019 is the shift from t/c to dx
# find and print t/c
x0 = foil[N/2].xc
c = foil[0].xc-x0
t = 2.*numpy.max([p.yc for p in foil])
#print "t/c = "+"%.3f"%(t/c)
# solve potential flow and boundary layer evolution
solve_gamma_kutta(foil,alpha)
top,bottom = solve_plot_boundary_layers(foil,alpha)
#Return point of seperation
return (top.x_sep-x0)/c, (bottom.x_sep-x0)/c
degrees = numpy.linspace(0, 10, 11) #For the graph axis
In [29]:
alpha = numpy.linspace(0, 10./180*numpy.pi, 11)
In [30]:
def split_panels(panels):
# positive velocity defines `top` BL
top = [p for p in panels if p.gamma<=0]
# negative defines the `bottom`
bottom = [p for p in panels if p.gamma>=0]
# reverse array so panel[0] is stagnation
bottom = bottom[::-1]
return top,bottom
In [31]:
sep_point = [] #Create empty list as it's easy.
for a in alpha:
sep_point.append(predict_jukowski_separation(0.15, a))
sep_point = numpy.array(sep_point) # Turn sep_point from list into an array
In [32]:
pyplot.figure(figsize=(10,7))
pyplot.ylabel(r'$\frac{x}{c}$', fontsize=24)
pyplot.xlabel(r'Angle of Attack $^o$', fontsize=16)
pyplot.plot(degrees, sep_point[:,0], label=r'Top Separation Point')
pyplot.plot(degrees, sep_point[:,1], label=r'Bottom Separation Point')
pyplot.legend(loc='right')
pyplot.show()
As we can see in the picture, when the angle of attack increases the top separation point move to the leading edge, however, the bottom separation point move to the trailing edge.
The geometry is the 15% Jukowski foil with camber, quantified by the height of the mean line at $x=0$, ie $\bar y = \frac 12 (y[N/4]+y[3N/4])$.
Compute the lift force coefficient $C_L$ and boundary layer separation point locations for $\bar y/c = 0.02,0.04,0.08$ at $\alpha=0$ using $N=128$ panels, and compare to the symmetric foil solution.
In [33]:
def camber(t_c, alpha, dy, N=128):
#Function to find empirical relationship between dy and camber.
dx=t_c-0.019
foil = make_jukowski(N, dx, dy)
#plot_flow(foil, alpha) #Can be commented in to inspect the shape of the foil
a = int(N/4)
b = int((3*N)/4)
y_bar = 0.5*(foil[a].yc + foil[b].yc)
return y_bar
trial_dy = numpy.linspace(0, -1, 30) #Range of experimental values for dy
y_bar = [] #Create empty list.
for i in trial_dy:
y_bar.append(camber(0.15, 0, i, 128))
In [34]:
#Investigating the relationship between dy and y_bar
def camber_emp(trial_dy, y_bar):
from sklearn.linear_model import LinearRegression
#This is a module developed for machine learning. It is useful for this sort of analysis.
X = trial_dy[:, numpy.newaxis]
regr = LinearRegression()
regr.fit(X, y_bar)
fig = pyplot.figure(figsize=(8,6))
ax = fig.add_subplot(1,1,1)
pyplot.ylabel(r"$\bar y$ ",fontsize=16)
pyplot.xlabel("dy ",fontsize=16)
pyplot.plot(trial_dy, y_bar)
pyplot.gca().set_xlim(left=-1)
pyplot.gca().set_ylim(bottom=0)
print("Multiply the dy value by %.5f to obtain the ybar/c value for the foil geometry chosen." % regr.coef_)
pyplot.show()
return regr.coef_
coef = camber_emp(trial_dy, y_bar)
we can see the realtionship between dy and y bar is linear correlation.
In [35]:
def jukowski_CL(alpha,t_c=0.15+0.019): return 2.*numpy.pi*(1+4/3/numpy.sqrt(3)*t_c)*numpy.sin(alpha)
def ybar_c(y_bar,c):
ybar_c=y_bar/c
return ybar_c
def predict_jukowski_separation_camber(t_c, ybar_c, alpha=0,N=128):
#Function from separation prediction and altered to account for camber.
# set dx to gets the correct t/c
dx = t_c -0.019
dy = ybar_c/(-0.84675)
foil = make_jukowski(N, dx, dy) #t_c-0.019 is the shift from t/c to dx
x0 = foil[N/2].xc
c = foil[0].xc-x0
t = 2.*numpy.max([p.yc for p in foil])
# solve potential flow and boundary layer evolution
solve_gamma_kutta(foil,alpha)
top,bottom = solve_plot_boundary_layers(foil,alpha)
#Return point of separation
return (top.x_sep-x0)/c, (bottom.x_sep-x0)/c
def calc_CL_camb(alpha, N, t_c, ybar_c):
#Function from earlier in this report altered to account for camber.
dx = t_c - 0.019
dy = ybar_c/(-0.84675)
foil = make_jukowski(N, dx, dy)
solve_gamma_kutta(foil, alpha)
return lift(foil)
fill_val = jukowski_CL(0) #The lift coefficient for the symmetric foil.
analytic = numpy.full_like(ybar_c, fill_val) #Python will get sad unless it gets an array of analytical C_L to plot.
analytic
Out[35]:
In [36]:
ybar_c = [0.02, 0.04, 0.08] #List of cambers.
c_l = []
for i in ybar_c:
c_l.append(calc_CL_camb(0, 128, 0.15, i))
sep_points_camber = []
for i in ybar_c:
sep_points_camber.append(predict_jukowski_separation_camber(0.15, i, 0, 128))
sep_points_camber = numpy.array(sep_points_camber)
In [ ]:
In [ ]:
In [ ]: